1 Setup

1.1 Load functions and set directories

Set knitting options

# global knitting options for automatic saving of all plots as .png and .pdf. Also sets cache directory.
knitr::opts_chunk$set(
  dev = c("png", "pdf"), fig.keep = "all",
  dev.args = list(pdf = list(encoding = "WinAnsi", useDingbats = FALSE)),
  fig.path = file.path("fig_output/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input()))),
  cache.path = file.path("cache/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input())))
)

Load packages

library(dada2); packageVersion("dada2") # sequence filtering and sample inference with dada2
## [1] '1.16.0'
library(ShortRead); packageVersion("ShortRead") # package for loading, displaying, and handling sequences
## [1] '1.46.0'
library(tidyverse); packageVersion("tidyverse") # handy piping and data operations
## [1] '1.3.0'
library(rjson); packageVersion("rjson") # read json output from figaro
## [1] '0.2.20'
# source all relevant scripting files
# paths.R contains paths to local data and programs
source(file.path("scripts", "paths.R"))

Set directories

# Set up names of sub directories to stay organized
preprocess_path <- file.path(project_path_OM16, "01_preprocess")
    demultiplex_path <- file.path(preprocess_path, "demultiplexed")
        demultiplex_stats_path <- file.path(demultiplex_path, "stats") # statistical summary output from idemp
        demultiplex_unassigned_path <- file.path(demultiplex_path, "unassigned")
        demultiplex_assigned_path <- file.path(demultiplex_path, "assigned")
    figaro_output_path <- file.path(preprocess_path, "Figaro_output")
filter_path <- file.path(project_path_OM16, "02_filter")

1.2 Organize and rename demultiplexing output files

(demux_output_raw_names <- sort(list.files(demultiplex_path, full.names = FALSE)))
## [1] "assigned"   "stats"      "unassigned"

Now we read in the names of the demultiplexed and assigned fastq files.

fnFs_raw_names <- sort(list.files(demultiplex_assigned_path, pattern = "R1", full.names = TRUE))
fnRs_raw_names <- sort(list.files(demultiplex_assigned_path, pattern = "R2", full.names = TRUE))

Extract sample names for later use

# Extract sample names
sample_names_F <- fnFs_raw_names %>% basename() %>% str_remove("Undetermined_S0_L001_R1_001.fastq.gz_") %>% str_remove(".fastq.gz")

sample_names_R <- fnRs_raw_names %>% basename() %>% str_remove("Undetermined_S0_L001_R2_001.fastq.gz_") %>% str_remove(".fastq.gz")

# Check that forward and reverse file names match
if(!identical(sample_names_F, sample_names_R)) stop("Forward and reverse file names do not match.")

# print first few names
sample_names_F %>% head()
## [1] "NSHQ10_2016" "NSHQ14_2016" "WAB103_2016" "WAB104_2016" "WAB105_2016"
## [6] "WAB12_2016"
# Repeat the above with reverse reads
R1_names_incremented <- vector("character", length(sample_names_F))
for (i in seq_along(R1_names_incremented )) {
  R1_names_incremented[[i]] <- paste0(sample_names_F[[i]], "_S", i, "_L001_R1_001.fastq.gz")
}

# Now, rename files
file.rename(from = fnFs_raw_names, to = file.path(demultiplex_assigned_path, R1_names_incremented))
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
fnFs <- sort(list.files(demultiplex_assigned_path, pattern = "R1", full.names = TRUE)) # re-assign full names of forward reads, now with Illumina naming convention for use with Figaro.

# Print first few lines of fnFs
fnFs %>% head()
## [1] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/NSHQ10_2016_S1_L001_R1_001.fastq.gz"
## [2] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/NSHQ14_2016_S2_L001_R1_001.fastq.gz"
## [3] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/WAB103_2016_S3_L001_R1_001.fastq.gz"
## [4] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/WAB104_2016_S4_L001_R1_001.fastq.gz"
## [5] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/WAB105_2016_S5_L001_R1_001.fastq.gz"
## [6] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned/WAB12_2016_S6_L001_R1_001.fastq.gz"
# Repeat the above with reverse reads
R2_names_incremented <- vector("character", length(sample_names_R))
for (i in seq_along(R2_names_incremented )) {
  R2_names_incremented[[i]] <- paste0(sample_names_R[[i]], "_S", i, "_L001_R2_001.fastq.gz")
}

file.rename(from = fnRs_raw_names, to = file.path(demultiplex_assigned_path, R2_names_incremented))
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
fnRs <- sort(list.files(demultiplex_assigned_path, pattern = "R2", full.names = TRUE))

Count and display reads for input

# count forward reads
input_read_count_F <- countLines(demultiplex_assigned_path, pattern="_R1_001.fastq")/4

# make forward read counts into tbl with sample name
input_read_count_F_tbl <- input_read_count_F %>% as_tibble(rownames = "file_name") %>% rename(reads_F = value) %>% mutate(sample_id = file_name %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz"))

# repeat for reverse

input_read_count_R <- countLines(demultiplex_assigned_path, pattern="_R2_001.fastq")/4

input_read_count_R_tbl <- input_read_count_R %>% as_tibble(rownames = "file_name") %>% rename(reads_R = value) %>% mutate(sample_id = file_name %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R2_001.fastq.gz"))


# Join forward and reverse together into one tbl
input_read_count_FR <- input_read_count_F_tbl %>% select(-file_name) %>% left_join(input_read_count_R_tbl %>% select(-file_name), by = "sample_id") %>% select(sample_id, everything())

# print. reads counts are same in forward and reverse, which is what we want.
input_read_count_FR

Plot read counts

plot_read_count_input <- input_read_count_FR %>% ggplot(aes(
  x = fct_reorder(sample_id, desc(reads_F)),
  y = reads_F
)) +
  geom_bar(stat = "identity") +
  scale_y_continuous(name = "Forward reads") +
  scale_x_discrete(name = "Sample ID") +
  theme_classic(base_size = 9)+
  theme(
    axis.text.x = element_text(angle = 90, vjust = .5, hjust = 1),
    legend.position = "bottom"
  )

plot_read_count_input

2 Check primer sequences

The 515 (forward) (Caporaso) and 806R (reverse) (Caporaso) primers (EarthMicrobiome) were used to amplify this dataset. We record the DNA sequences for those primers.

However, the sequencing strategy used here is not supposed to sequence the primers.

FWD <- "GTGCCAGCMGCCGCGGTAA" # 515F (Caporaso)
REV <- "GGACTACHVGGGTWTCTAAT" # 806R (Caporaso)

(FWD_length <- str_length(FWD))
## [1] 19
(REV_length <- str_length(REV))
## [1] 20

That said, let’s check for any primer sequence matches in the sequences.

allOrients <- function(primer) {
    # Create all orientations of the input sequence
    require(Biostrings)
    dna <- DNAString(primer)  # The Biostrings works w/ DNAString objects rather than character vectors
    orients <- c(Forward = dna, Complement = complement(dna), Reverse = reverse(dna), 
        RevComp = reverseComplement(dna))
    return(sapply(orients, toString))  # Convert back to character vector
}
FWD.orients <- allOrients(FWD)
REV.orients <- allOrients(REV)
FWD.orients
##               Forward            Complement               Reverse 
## "GTGCCAGCMGCCGCGGTAA" "CACGGTCGKCGGCGCCATT" "AATGGCGCCGMCGACCGTG" 
##               RevComp 
## "TTACCGCGGCKGCTGGCAC"

We are now ready to count the number of times the primers appear in the forward and reverse read, while considering all possible primer orientations. Identifying and counting the primers on one set of paired end FASTQ files is sufficient, assuming all the files were created using the same library preparation, so we’ll just process the first sample.

# write function to check for primer sequences
primerHits <- function(primer, fn) { # Counts number of reads in which the primer is found
    nhits <- vcountPattern(primer, sread(readFastq(fn)), fixed = FALSE)
    return(sum(nhits > 0))
}
# print file name
print(paste("file_name: ", basename(fnFs[[1]])))
## [1] "file_name:  NSHQ10_2016_S1_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs[[1]])/4))
## [1] "total reads: 26412"
# check for primers
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs[[1]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnRs[[1]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs[[1]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnRs[[1]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0
# print file name
print(paste("file_name: ", basename(fnFs[[2]])))
## [1] "file_name:  NSHQ14_2016_S2_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs[[2]])/4))
## [1] "total reads: 20360"
# check for primers
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs[[2]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnRs[[2]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs[[2]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnRs[[2]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       1
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0

There were essentially no primer sequences detected. We will check if primer sequences are still found later after we truncate lower quality reads. So now let’s inspect read quality profiles.

3 Inspect read quality profiles

It’s important to get a feel for the quality of the data that we are using. To do this, we will plot the quality of some of the samples.

From the dada2 tutorial: >In gray-scale is a heat map of the frequency of each quality score at each base position. The median quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position (this is more useful for other sequencing technologies, as Illumina reads are typically all the same lenghth, hence the flat red line).

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs) <= 20) {
  plotQualityProfile(fnFs)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs)) # grab 20 random samples to plot
  plotQualityProfile(fnFs[rand_samples])
}

# Reverse reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnRs) <= 20) {
  plotQualityProfile(fnRs)
} else {
  rand_samples <- sample(size = 20, 1:length(fnRs)) # grab 20 random samples to plot
  plotQualityProfile(fnRs[rand_samples])
}

4 Filter parameter optimization with Figaro

Create a directory for Figaro output

# Create directory to hold the output from figaro
if (!dir.exists(figaro_output_path)) dir.create(figaro_output_path)

I couldn’t get figaro (v1.1.1) to run by system2() in R due to a shell permission denied error, but it worked on the command line. This chunk generates the exact text I entered in the command line after changing to the figaro directory. I then just read in the figaro output JSON and continue from there.

# text to change directory in terminal to figaro 
paste("cd", figaro_path) %>% str_remove("/figaro.py") %>% cat()
## cd /Users/melo.d/opt/figaro
# primers were not sequenced
# sequences are 1X160 bp
# Tried 240, lost nearly all reads at merge step.
# Try again with 260
figaro_amplicon_length <- 260

figaro_min_overlap <- 20

figaro_output_file_name <- "Figaro_params_OM16.json"

# using 1 for primer lengths because even though primers were not sequenced, the current version of figaro does not allow 0 as a primer length.

flags_figaro <- paste("--ampliconLength", figaro_amplicon_length, "--forwardPrimerLength", 1, "--reversePrimerLength", 1, "--outputFileName", figaro_output_file_name, "--inputDirectory", demultiplex_assigned_path, "--outputDirectory", figaro_output_path, "--minimumOverlap", figaro_min_overlap)

paste("python3", "figaro.py", flags_figaro) %>% cat()
## python3 figaro.py --ampliconLength 260 --forwardPrimerLength 1 --reversePrimerLength 1 --outputFileName Figaro_params_OM16.json --inputDirectory /Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/demultiplexed/assigned --outputDirectory /Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/01_preprocess/Figaro_output --minimumOverlap 20

Read in Figaro JSON output, top 5 suggested trim parameters.

# Give the input file name to the function.
figaro_result <- fromJSON(file = file.path(figaro_output_path, figaro_output_file_name))

figaro_result[1:5] %>% paste(collapse = '\n') %>% cat()
## list(trimPosition = c(132, 150), maxExpectedError = c(1, 1), readRetentionPercent = 93.14, score = 93.1438541871584)
## list(trimPosition = c(133, 149), maxExpectedError = c(1, 1), readRetentionPercent = 93.13, score = 93.1327351882104)
## list(trimPosition = c(134, 148), maxExpectedError = c(1, 1), readRetentionPercent = 93.13, score = 93.1301692653763)
## list(trimPosition = c(135, 147), maxExpectedError = c(1, 1), readRetentionPercent = 93.03, score = 93.0318088900673)
## list(trimPosition = c(136, 146), maxExpectedError = c(1, 1), readRetentionPercent = 92.96, score = 92.9599630507112)
# Save and print Figaro-suggested parameters for dada2::filterAndTrim
(truncLen_best_figaro <- figaro_result[[1]]$trimPosition)
## [1] 132 150
(maxEE_best_figaro <- figaro_result[[1]]$maxExpectedError)
## [1] 1 1

5 Filter and trim for quality

5.1 Organize file paths and directories

# set file path for quality-filtered read
filter_path <- file.path(project_path_OM16, "02_filter") 

# Put filtered reads into separate sub-directories for big data workflow

if (!dir.exists(filter_path)) dir.create(filter_path)
    subF_path <- file.path(filter_path, "preprocessed_F") 
    subR_path <- file.path(filter_path, "preprocessed_R") 
if (!dir.exists(subF_path)) dir.create(subF_path)
if (!dir.exists(subR_path)) dir.create(subR_path)
    
# Move demultiplexed and assigned R1 and R2 to separate forward/reverse sub-directories

fnFs.Q <- file.path(subF_path,  basename(list.files(demultiplex_assigned_path , pattern="R1", full.names = TRUE))) 
fnRs.Q <- file.path(subR_path,  basename(list.files(demultiplex_assigned_path , pattern="R2", full.names = TRUE)))

file.rename(from = list.files(demultiplex_assigned_path , pattern="R1", full.names = TRUE), to = fnFs.Q)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
file.rename(from = list.files(demultiplex_assigned_path , pattern="R2", full.names = TRUE), to = fnRs.Q)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# File parsing; create file names and make sure that forward and reverse files match
filtpathF <- file.path(subF_path, "filtered") # files go into preprocessed_F/filtered/
filtpathR <- file.path(subR_path, "filtered") # ...
fastqFs <- sort(list.files(subF_path, pattern="fastq.gz"))
fastqRs <- sort(list.files(subR_path, pattern="fastq.gz"))
if(length(fastqFs) != length(fastqRs)) stop("Forward and reverse files do not match.")

Before chosing sequence variants, we want to trim reads where their quality scores begin to drop (the truncLen and truncQ values) and remove any low-quality reads that are left over after we have finished trimming (the maxEE value).

5.2 Filter and trim

# No argument pasted to trimLeft because the primers were not sequenced
filt_out <- filterAndTrim(fwd = file.path(subF_path, fastqFs), filt = file.path(filtpathF, fastqFs), rev = file.path(subR_path, fastqRs), filt.rev = file.path(filtpathR, fastqRs), truncLen = truncLen_best_figaro, maxEE = maxEE_best_figaro, truncQ = 2, maxN = 0, rm.phix = TRUE, compress = TRUE, verbose = TRUE, multithread = TRUE)

# look at how many reads were kept
head(filt_out)
##                                     reads.in reads.out
## NSHQ10_2016_S1_L001_R1_001.fastq.gz    26412     24508
## NSHQ14_2016_S2_L001_R1_001.fastq.gz    20360     19195
## WAB103_2016_S3_L001_R1_001.fastq.gz    21379     19862
## WAB104_2016_S4_L001_R1_001.fastq.gz    22421     20660
## WAB105_2016_S5_L001_R1_001.fastq.gz    28289     26099
## WAB12_2016_S6_L001_R1_001.fastq.gz     26308     24403
# summary of samples in filt_out by percentage
filt_out %>% 
  data.frame() %>% 
  mutate(Samples = rownames(.),
         percent_kept = 100*(reads.out/reads.in)) %>%
  select(Samples, everything()) %>%
  summarise(min_remaining = paste0(round(min(percent_kept), 2), "%"), 
            median_remaining = paste0(round(median(percent_kept), 2), "%"),
            mean_remaining = paste0(round(mean(percent_kept), 2), "%"), 
            max_remaining = paste0(round(max(percent_kept), 2), "%"))

5.3 Check if primers were removed by filterAndTrim

Generate sorted file name vector for filtered files

fnFs_filt <- sort(list.files(filtpathF, full.names = TRUE))

# print first 6 lines
fnFs_filt %>% head()
## [1] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/NSHQ10_2016_S1_L001_R1_001.fastq.gz"
## [2] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/NSHQ14_2016_S2_L001_R1_001.fastq.gz"
## [3] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/WAB103_2016_S3_L001_R1_001.fastq.gz"
## [4] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/WAB104_2016_S4_L001_R1_001.fastq.gz"
## [5] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/WAB105_2016_S5_L001_R1_001.fastq.gz"
## [6] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2016/02_filter/preprocessed_F/filtered/WAB12_2016_S6_L001_R1_001.fastq.gz"
fnRs_filt <- sort(list.files(filtpathR, full.names = TRUE))
# print file name
print(paste("file_name: ", basename(fnFs_filt[[1]])))
## [1] "file_name:  NSHQ10_2016_S1_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs_filt[[1]])/4))
## [1] "total reads: 24508"
# tally primer hits in file
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[1]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[1]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[1]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[1]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0
# print file name
print(paste("file_name: ", basename(fnFs_filt[[2]])))
## [1] "file_name:  NSHQ14_2016_S2_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs_filt[[2]])/4))
## [1] "total reads: 19195"
# tally primer hits in file
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[2]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[2]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[2]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[2]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0

Primer sequences appear to be absent. Thus, additional filtering of primers is probably not needed, and won’t be performed here.

5.4 Plot the quality of the filtered and cut fastq files.

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs_filt) <= 20) {
  plotQualityProfile(fnFs_filt)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs_filt)) # grab 20 random samples to plot
  plotQualityProfile(fnFs_filt[rand_samples])
}

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs_filt) <= 20) {
  plotQualityProfile(fnFs_filt)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs_filt)) # grab 20 random samples to plot
  plotQualityProfile(fnFs_filt[rand_samples])
}

6 Infer sequence variants

In this part of the pipeline, dada2 will learn to distinguish error from biological differences using a subset of our data as a training set. After it understands the error rates, we will reduce the size of the dataset by combining all identical sequence reads into “unique sequences”. Then, using the dereplicated data and error rates, dada2 will infer the sequence variants in our data. Finally, we will merge the corresponding forward and reverse reads to create a list of the fully denoised sequences and create a sequence table from the result.

6.1 Learn the error rates

set.seed(100) # set seed to ensure that randomized steps are reproducible

# Learn forward error rates
errF <- learnErrors(fnFs_filt, nbases=1e8, multithread=TRUE)
## 28671456 total bases in 217208 reads from 10 samples will be used for learning the error rates.
# Learn reverse error rates
errR <- learnErrors(fnRs_filt, nbases=1e8, multithread=TRUE)
## 32581200 total bases in 217208 reads from 10 samples will be used for learning the error rates.

6.1.1 Plot Error Rates

The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. Here, we want to see estimated error rates (black line) that are a good fit to the observed rates (points), and that the error rates drop with increased quality as expected. If these criteria aren’t observed, then it may be a good idea to increase the nbases parameter. This allows the machine learning algorithm to train on a larger portion of your data and may help improve the fit.

plotErrors(errF, nominalQ=TRUE)
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

plotErrors(errR, nominalQ=TRUE)
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

6.2 Dereplication, sequence inference, and merging of paired-end reads

Extract sample names and add them to vector of filtered file names

# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample_names_F_filt <- fnFs_filt %>% basename() %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz")
sample_names_R_filt <- fnRs_filt %>% basename() %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R2_001.fastq.gz")

# Double check
if(!identical(sample_names_F_filt, sample_names_R_filt)) stop("Forward and reverse files do not match.")

names(fnFs_filt) <- sample_names_F_filt
names(fnRs_filt) <- sample_names_R_filt
# make lists to hold the loop output
mergers <- vector("list", length(sample_names_F))
names(mergers) <- sample_names_F
ddF <- vector("list", length(sample_names_F))
names(ddF) <- sample_names_F
ddR <- vector("list", length(sample_names_F))
names(ddR) <- sample_names_F

# For each sample, get a list of merged and denoised sequences
for(sam in sample_names_F) {
    cat("Processing:", sam, "\n")
    # Dereplicate forward reads
    derepF <- derepFastq(fnFs_filt[[sam]])
    # Infer sequences for forward reads
    ddF[[sam]] <- dada(derepF, err=errF, multithread=TRUE)
    # Dereplicate reverse reads
    derepR <- derepFastq(fnRs_filt[[sam]])
    # Infer sequences for reverse reads
    ddR[[sam]] <- dada(derepR, err=errR, multithread=TRUE)
    # Merge reads together
    merger <- mergePairs(ddF[[sam]], derepF, ddR[[sam]], derepR)
    mergers[[sam]] <- merger
}
## Processing: NSHQ10_2016 
## Sample 1 - 24508 reads in 4174 unique sequences.
## Sample 1 - 24508 reads in 4814 unique sequences.
## Processing: NSHQ14_2016 
## Sample 1 - 19195 reads in 2999 unique sequences.
## Sample 1 - 19195 reads in 3358 unique sequences.
## Processing: WAB103_2016 
## Sample 1 - 19862 reads in 3920 unique sequences.
## Sample 1 - 19862 reads in 4053 unique sequences.
## Processing: WAB104_2016 
## Sample 1 - 20660 reads in 5493 unique sequences.
## Sample 1 - 20660 reads in 5597 unique sequences.
## Processing: WAB105_2016 
## Sample 1 - 26099 reads in 6489 unique sequences.
## Sample 1 - 26099 reads in 7005 unique sequences.
## Processing: WAB12_2016 
## Sample 1 - 24403 reads in 3971 unique sequences.
## Sample 1 - 24403 reads in 4015 unique sequences.
## Processing: WAB188_2016 
## Sample 1 - 22989 reads in 5194 unique sequences.
## Sample 1 - 22989 reads in 5381 unique sequences.
## Processing: WAB55_2016 
## Sample 1 - 16081 reads in 3526 unique sequences.
## Sample 1 - 16081 reads in 3544 unique sequences.
## Processing: WAB56_2016 
## Sample 1 - 21403 reads in 3412 unique sequences.
## Sample 1 - 21403 reads in 3771 unique sequences.
## Processing: WAB71_2016 
## Sample 1 - 22008 reads in 3511 unique sequences.
## Sample 1 - 22008 reads in 3835 unique sequences.
rm(derepF); rm(derepR)

6.3 Construct sequence table

seqtab <- makeSequenceTable(mergers)
write_rds(seqtab, path = format(Sys.Date(), "data_output/seqtab_OM16_processed_%Y%m%d.rds"), compress = "gz")

7 Remove Chimeras

Although dada2 has searched for insertion/deletion errors and substitutions, there may still be chimeric sequences in our dataset (sequences that are derived from forward and reverse sequences from two different organisms becoming fused together during PCR and/or sequencing). To identify chimeras, we will search for rare sequence variants that can be reconstructed by combining left-hand and right-hand segments from two more abundant “parent” sequences.

# Remove chimeras
seqtab_nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE)
# Print percentage of our sequences that were not chimeric.
100*sum(seqtab_nochim)/sum(seqtab)
## [1] 98.35169
write_rds(seqtab_nochim, path = format(Sys.Date(), "data_output/seqtab_nochim_OM16_processed_%Y%m%d.rds"), compress = "gz")

8 Track reads through pipeline

getN <- function(x) sum(getUniques(x)) # function to grab sequence counts from output objects
# tracking reads by counts
filt_out_track <- filt_out %>%
  data.frame() %>%
  mutate(Sample = rownames(.) %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz")) %>%
  rename(input = reads.in, filtered = reads.out)
rownames(filt_out_track) <- filt_out_track$Sample

ddF_track <- data.frame(denoisedF = sapply(ddF[sample_names_F_filt], getN)) %>%
  mutate(Sample = row.names(.))
ddR_track <- data.frame(denoisedR = sapply(ddR[sample_names_F_filt], getN)) %>%
  mutate(Sample = row.names(.))
merge_track <- data.frame(merged = sapply(mergers, getN)) %>%
  mutate(Sample = row.names(.))
chim_track <- data.frame(nonchim = rowSums(seqtab_nochim)) %>%
  mutate(Sample = row.names(.))

track <- left_join(filt_out_track, ddF_track, by = "Sample") %>%
  # left_join(ddF_track, by = "Sample") %>%
  left_join(ddR_track, by = "Sample") %>%
  left_join(merge_track, by = "Sample") %>%
  left_join(chim_track, by = "Sample") %>%
  replace(., is.na(.), 0) %>%
  select(Sample, everything())
row.names(track) <- track$Sample
head(track)
# tracking reads by percentage
track_pct <- track %>% 
  data.frame() %>%
  mutate(Sample = rownames(.),
         filtered_pct = ifelse(filtered == 0, 0, 100 * (filtered/input)),
          # cut_pct = ifelse(cut == 0, 0, 100 * (cut/filtered)),
         denoisedF_pct = ifelse(denoisedF == 0, 0, 100 * (denoisedF/filtered)),
         denoisedR_pct = ifelse(denoisedR == 0, 0, 100 * (denoisedR/filtered)),
         merged_pct = ifelse(merged == 0, 0, 100 * merged/((denoisedF + denoisedR)/2)),
         nonchim_pct = ifelse(nonchim == 0, 0, 100 * (nonchim/merged)),
         total_pct = ifelse(nonchim == 0, 0, 100 * nonchim/input)) %>%
  select(Sample, ends_with("_pct"))

# summary stats of tracked reads averaged across samples
track_pct_avg <- track_pct %>% summarize_at(vars(ends_with("_pct")), 
                                            list(avg = mean))
head(track_pct_avg)
track_pct_med <- track_pct %>% summarize_at(vars(ends_with("_pct")), 
                                            list(avg = stats::median))
head(track_pct_avg)
head(track_pct_med)
# Plotting each sample's reads through the pipeline
track_plot <- track %>% 
  pivot_longer(-Sample, names_to = "Step", values_to = "Reads") %>% 
  mutate(Step = factor(Step, 
                       levels = c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim"))) %>%
  ggplot(aes(x = Step, y = Reads)) +
  geom_line(aes(group = Sample), alpha = 0.2) +
  geom_point(alpha = 0.5, position = position_jitter(width = 0)) + 
  stat_summary(fun = median, geom = "line", group = 1, color = "steelblue", size = 1, alpha = 0.5) +
  stat_summary(fun = median, geom = "point", group = 1, color = "steelblue", size = 2, alpha = 0.5) +
  stat_summary(fun.data = median_hilow, fun.args = list(conf.int = 0.5), 
               geom = "ribbon", group = 1, fill = "steelblue", alpha = 0.2) +
  geom_label(data = t(track_pct_avg[1:5]) %>% data.frame() %>% 
               rename(Percent = 1) %>%
               mutate(Step = c("filtered", "denoisedF", "denoisedR", "merged", "nonchim"),
                      Percent = paste(round(Percent, 2), "%")),
             aes(label = Percent), y = 1.1 * max(track[,2])) +
  geom_label(data = track_pct_avg[6] %>% data.frame() %>%
               rename(total = 1),
             aes(label = paste("Total\nRemaining:\n", round(track_pct_avg[1,6], 2), "%")), 
             y = mean(track[,6]), x = 6.6) +
  expand_limits(y = 1.1 * max(track[,2]), x = 7.2) +
  theme_classic()

track_plot